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CONTRACTOR REPORT 


COMPUTATION OF TURBULENT BOUNDARY LAYER FLOWS WITH 
AN ALGEBRAIC STRESS TURBULENCE MODEL 

INTRODUCTION 


During the last decades, significant progress has been made in mathematically 
modeling the turbulent flows. Memorable events in the effort, among many others, 
would be the two Stanford Conferences on turbulent flows [1,2]. 

Development of turbulence models, and accordingly computation of turbulent 
flows too, were made, entirely to certain extents, using one or another form of the 
finite difference computational method. Only very recently, a few publications on 
finite element computation of turbulent flows begin to appear. 

A list of publications on finite element computation of turbulent flows can be 
found in Reference 3. To summarize, finite element computation of turbulent flows 
resulted in, more or less, dissatisfaction in most of the cases. The turbulent kinetic 
energy or dissipation rate becomes negative during iterative solution procedure so that 
a dedicated computational procedure to overcome this difficulty has also been proposed 
[ 4] . The hardly obtained computational results were less satisfactory than that of the 
finite difference computation of the same turbulent flow using the same turbulence 
model in some cases [ 4] . Considering the success of the finite element method in 
other areas of engineering computations, such as the structural mechanics, laminar 
flow problems, and chemical process modelings to mention a few, it is quite unusual 
that the finite element method can hardly be used for computation of turbulent flows 
especially when k-e type turbulence models are used. But it should be mentioned 
here that there already exist commercial finite element flow analysis codes using the 
mixing length turbulence models [ 5] . 

A mathematically rigorous finite element turbulent boundary layer flow analysis 
code has been constructed, which is based on the finite element laminar boundary 
layer flow analysis code, to assess and/or validate turbulence models to be used for 
general finite element computation of turbulent flows [3]. 

It was found that the flow field variables (velocity, turbulent kinetic energy, 
and dissipation rate) can hardly make smooth transition from the near wall high 
turbulence region to the free stream low turbulence region in most of the finite 
element computations. This difficulty was able to be removed by employing a type 
of algebraic stress turbulence model. The remaining propositions included in the 
present turbulence model contribute to improving the computational results for turbu- 
lent boundary layer flows. Use of the present turbulence model in finite difference 
computation of turbulent boundary layer flows, as well as elliptic flows, improved the 
computational results compared with those obtained by using the standard k-e 
turbulence models [ 7] . 

It would be worthwhile to mention that most of the k-e turbulence models were 
able to run smoothly in the present finite element turbulent boundary layer flow code 
(BLFLOW) [3] simply by including the algebraic stress form of the eddy viscosity into 
those turbulence models. In this regard, the observed differences between the two 
numerical methods are also discussed in this report. 


TURBULENCE MODEL 


Eddy Viscosity Equation 

The eddy viscosity expression can be obtained by contracting the Reynolds 
stress turbulence model [8] using an algebraic Reynolds stress transport assumption 
[9], In most of the algebraic stress turbulence models, the eddy viscosity is given 
as [9,10] : 
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where v, is the turbulent eddy viscosity; k is the turbulent kinetic energy; e is the 
^ 2 

dissipation rate; P, P = v^(9u/3y) , is the production rate; u is the time averaged 

flow direction velocity; y is the transverse coordinate; c ^ = 3.7 + 0.7 tanh(P/e); 

and c „ = 0.32. ^ 
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The functional form of c^^ was used to bring down c^ to 0.09 for equilibrium 

state (production rate is equal to the dissipation rate of turbulent kinetic energy) . 
The c^ curve is shown in Figure 1 together with that of Launder [ 10] . A different 


c^ curve can also be found in Rodi [9]. Experimentally deduced c^ curves [9] are 
also shown in Figure 1. It can be seen that the values of at the center region of 


jets are close to those at equilibrium state for both jets, respectively. The present 
c^ function cannot model this phenomenon. The authors are not aware of any turbu- 
lence model which can achieve this phenomenon. This defect deteriorates part of the 
computational results, which is discussed in the following section. Nevertheless, it 
was found that the c^ function was definitely required for the present finite element 


computation of turbulent boundary layer flows to ensure smooth transition of flow 
field variables (u, k, e) from the near wall high turbulence region to the free stream 
low turbulence region. For multi- dimensional elliptic flows, the potential core region 
constitutes the low turbulence region, and the c function is required for the same 
reason as before. ^ 


It can be seen in Reference 7 that the c^ function is not required for the 

finite difference computation of turbulent flows, and that the same turbulence model 
with constant c^ (=0.09) performed much better than the present case. 


Dissipation Rate Equation 

It has long been recognized that the standard dissipation rate equation is lack 
of variable energy transfer capability which will transfer more turbulent kinetic energy 
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into dissipation range when production is high. The same line of efforts can be 
found in Pope [11], Hanjalic and Launder [12], and in Hanjalic et al., [13] among 
many others. 

In the present study, it is proposed to let part of the dissipation rate (e) 
evolve according to a time- scale related to the production itself. The proposed dissi- 
pation rate equation is given as: 
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where D/Dt represents material derivative; v is the kinematic viscosity, r is an integer 
constant such that r = 0 for two-dimensional flows and r = 1 for axisymmetric flows; 
and a is a constant. In equation (3), k/e can be considered to be a time scale 
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related to the dissipation rate; and k/P, to the production rate. After simplification, 
equation (3) becomes: 
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As usual, (Cj^+Cg) may be determined based on the generation rate of dissipation 

rate as given in Harris et al. [14] and C 2 may be determined according to the decay 

of grid turbulence experiment of Harlow and Nakayama [15]. The constants used in 
the present study are: Cj^ = 1.15, C 2 = 1.90, and Cg = 0.25. 


Diffusion Coefficients for Turbulence Quantities 
Transport equation for turbulent kinetic energy is given as: 
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Many of the experimental data, such as the flat plate boundary layer flow due 
to Klebanoff [ 16] and turbulent jets exhausting into moving streams considered herein 
as example problems among many others, shows that the turbulent kinetic energy and 
its gradient are persistently big even after the velocity gradient has vanished near 
and in the free stream region. It is conjectured that the turbulent kinetic energy has 
been diffused outward from the high turbulence region due to the following reasons. 
Firstly, the velocity gradient has vanished and there is no production of turbulent 
kinetic energy in the region. Secondly, the boundary layer flows are spreading in 
the downstream direction, hence the turbulent kinetic energy could not have been 
convected from the upstream region. Thirdly, the decay rate of turbulence, as can 
be found in the grid turbulence, is so small that the same level of turbulent kinetic 
energy should exist throughout the free stream region if it has been convected from 
the high turbulence region. Consequently, the diffusivity of turbulent kinetic energy 
needs to be bigger than that of the momentum equation, and hence the diffusion 
coefficient, for turbulent kinetic energy was set to be smaller than unit to 

achieve the observed diffusion process of the turbulent kinetic energy. The proposed 

values for the model constants are a, =0.75 and a =1.05. The value of o 

k e e 

definitely needs to be bigger than for realizability requirement, i.e., the dissipa- 
tion rate should decay faster than the turbulent kinetic energy as the free stream 
region is approached. It was found that a, = 1.0 and a =1.3, which are most 

K S 

frequently used values in finite difference computations, were acceptable for the 
present finite element computation even though the diffusion process of turbulent 
kinetic energy near the free stream region was not so desirable. 
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TURBULENT FLOW EQUATIONS 


For the sake of clarity and completeness, the standard turbulent boundary layer 
equations are given below. 


^ + J: L 

3 X r 3y 

y 


(y"^ V) 


0 


( 6 ) 


3u IH - -1- _L 
3x ''^3y~yr 3y 


y^ (v + 


vt) 



dx 


(7) 


where u and v are time averaged velocities in flow direction and in transverse direc- 
tion, respectively, and the rest of the notations are the same as before. 


FINITE ELEMENT COMPUTATIONAL PROCEDURE 


Details of the present finite element computational procedure for turbulent 
boundary layer flows can be found in References 3 and 6. The coarse grid accuracy 
reported in Reference 6 was found to hold for the present turbulent boundary layer 
flow computations too [ 6] . Only the specific aspects of the present method and the 
fundamental differences between the finite element [3] and the finite difference 
method [7] as observed by the authors are described below. 

The transverse computational domain was extended outward several times of the 
boundary layer thicknesses at the initial line of the flow direction computational 
domain, and orthogonal grids were used everywhere in the flow domain in order to 
avoid even the slightest source of numerical uncertainty that can be caused by non- 
orthogonal grid transformation. This scheme is found to be especially helpful in 
assessing turbulence models to be used for multi- dimensional elliptic flow computations 
[ 3] . Since the low turbulence region almost always exists in elliptic flows in the form 
of potential core, it is important to conform smooth transition of flow field variables 
from the high turbulence region to low turbulence region. 

For linear elliptic boundary value problems, the convergence rate is given as 

[17]: 


Hello = Ch^^^ (8) 

where j |ej 1^, 1 le| |q ~ zeroth order error norm; is the com- 

putational domain; a is the exact solution; a^^ is the finite element solution; C is a 

constant dependent on problems; h is the mesh size; and k is the order of interpolat- 
ing polynomial. Equation (6) denotes that the convergence rate of the zeroth order 
error norm would be 3 if quadratic elements are used. For parabolic differential 
equations, the error norm is bounded by the discretization error in time domain and 
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that of the spatial domain [18]. For the boundary layer flows, the previous compu- 
tational experiment showed that the convergence rate is approximately 2.6 for fixed 
step-size (Ax) in flow direction [6]. This degeneration of convergence rate was 
caused by the discretization error in the flow -direction (time-like domain) as well as 
the feed-back effect of the non-linear time-like derivative, u(3u/9x). Nevertheless, 
considering that the error norm is defined in global sence, that the convergence rate 
is still higher than two, and that h is the real mesh size used in computations, it 
may be claimed that the quadratic finite element computation is more accurate than a 
second-order-accurate finite differencing scheme for which the convergence rate is 
obtained from the Taylor series expansion as the mesh size becomes infinitely small. 
This statement can be supported by many available computational results on heat con- 
duction problems and laminar flow calculations. Based on these facts, it can also be 
claimed that the numerical diffusion involved in the finite element computation would 
be less compared to a second-order- accurate finite differencing scheme. This can 
partly explain why the experimentally observed c^ -variation in the turbulent flow 

field needs to be incorporated into the finite element computation in order to achieve 
the realizability (i.e., the turbulent kinetic energy and the dissipation rate need to 
be positive) of turbulence quantities. 


COMPUTATIONAL RESULTS 


All of the computations were performed on the physical domain using physical 

dimensions; and the material property data used are the density, p, equal to 1.225 

3 - 4 

kg/m , and the molecular viscosity, y , equal to 0.17854 x 10 kg/m-sec. For all of 

the example cases, 45 unequally spaced quadratic elements were used to discretize the 

transverse domain, unless otherwise specified. The convergence criterion used is 

given as: 
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< 1 X 10 


1,N 
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where the superscript n denotes iteration level, the subscript j denotes each degree 
of freedom, and N denotes the total degree of freedom. 


1. Fully Developed Channel Flow 

Experimental data for the fully developed channel flow used herein can be found 
in Laufer [19]. Half of the channel width is equal to 0.0625 m, the center line mean 
velocity is equal to 7.07 m/sec, the Reynolds number based on these two parameters 

2 3 

is approximately 30,800, the pressure gradient is equal to -1.405 kg-m/sec -m , and 
the wall friction velocity, u^ is equal to 0.270 m/sec. 

The computational domain which extends from y = 0.005 m near the wall, which 
corresponds to y'*' = 100, to y = 0.0635 m at the center of the channel was discretized 
using 20 equally spaced quadratic elements. Near y"*" = 40 and below, high Reynolds 
number turbulence models may not yield good predictions unless these are modified by 
including low turbulence Reynolds number effects. 
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The Dirichlet boundary conditions for u and k at the near wall region were 
obtained from experimental data, whereas the boundary condition for e was obtained 
from the mixing length assumption given below. 
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where SL is the turbulence length scale, is the mixing length, k = 0.41, and c^^ = 
0.09. The near wall boundary conditions used are: u = 5.084 m/sec; k = 0.213 
m^/sec^; and e = 10.64 m^/sec^. 


A total of 84 iterations were used to achieve the convergence criterion given 
previously. The computational results for the velocity, the turbulent kinetic energy 
and the Reynolds stress are compared with experimental data in Figure 2. 



X: VELOCITY (u),0: KINETIC ENERGY Ik), a: SHEARING STRESS |-u'v') 
. : COMPUT. RESULT FOR u, k, AND -u'v' RESPECTIVELY. 


Figure 2. Fully developed channel flow. 
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2. Fully Developed Pipe Flow 

The fully developed pipe flow is another case of elliptic boundary layer flow. 
The experimental data used can be found in Laufer [ 20] . The diameter of the pipe is 
0.24688 m, the center line velocity is 30.48 m/sec, the Reynolds number based on 
these two parameters is 500,000, the pressure gradient in the flow direction is equal 

2 3 

to -23.05 kg-m/sec -m , and the wall friction velocity is equal to 1.078 m/sec. 

The computational domain extending from y = 0 m at the center of the pipe to 
y = 0.1192 m at the near wall region, which corresponds to y'*' = 300, were discretized 
using 20 equally spaced quadratic elements. The near wall boundary conditions 

2 2 

obtained from experimental data are given as: u = 20 m/sec, k = 4.93 m /sec , and 

2 3 

e = 1000 m /sec . It took 100 iterations to achieve the same prescribed convergence 
criterion. The computational results are compared with experimental data in Figure 3. 
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NOTATIONS ARE THE SAME AS IN FIGURE 1. 

Figure 3. Fully developed pipe flow. 
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3. Flat Plate Boundary Layer Flow Under Vanishing Pressure Gradient 

The Wiehardt and Tillman flat plate boundary layer flow [21] is considered 
below. Another case of equilibrium boundary layer flow with zero pressure gradient 
can be found in Klebanoff [ 16] . Both of the flat plate boundary layer flows are self- 
similar flows. Experimental wall shearing stresses are presented in Wiehardt and 
Tillmann [21], whereas self-similar Reynolds stresses are given in Klebanoff [16]. 

In the following computation, the initial condition data for velocity were obtained from 
Wiehardt and Tillmann [21]; and the turbulent kinetic energy from Klebanoff [16]. 

The free stream velocity is equal to 33 m/sec and the free stream turbulence 
level is 0.25 percent. The transverse computational domain extending from y = 0.001 
m at the near waU region, which corresponds to y'*' = 100, to y = 0.112 m at the free 
stream region was discretized using 45 unequally spaced quadratic elements. The 
computational domain in the flow direction, which extends from x = 0.937 m up to x = 
4.987 m was discretized using 600 line-steps. 

Development of the wall shearing stress along the flow direction is shown in 
Figure 4; the computed velocity profile, the kinetic energy profile, and the Reynolds 
stress are compared with experimental data in Figure 5. Approximately 14 iterations 
were required for each line- step. 



Figure 4. V^all shearing stress for the flat plate flow. 
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4. A Plane Jet Exhausting into a Moving Stream 

A plane jet exhausting into a moving stream is considered below, the experi- 
mental data of which can be found in Bradbury [22]. 

The transverse domain extends from center line of the jet, y = 0, to y = 0.1 m 
which is approximately 12.5 times of the half jet width (the half jet width is defined 
as a distance from the center line of jet to a location where the excess velocity is half 
of the center line excess velocity). The computational domain in flow direction starting 
from X = 0.095 m (x/d = 10, d = 0.009525 m is the jet exit width) to x = 0.65 m 
(x/d = 70) was discretized by 580 line- steps. 

Decay of the center line velocity and evolution of the half jet width along the 
flow direction compares favorably with experimental data as shown in Figure 6. The 
computed velocity profile, turbulent kinetic energy profile, and the Reynolds stress 
are compared with experimental data in Figure 7 at two flow direction locations. 

These profiles compare favorably with experimental data at the middle of the flow 
direction domain, whereas these profiles compare less favorably with experimental data 
at the end of the computational domain. These degenerated solutions may be due to 
less appropriate c^ function. But the global features such as the growth of the half 

jet width were found to be less sensitive to the c^ function. 

An average of 11 iterations were required to satisfy the convergence criterion 
for each line-step. 


o' 

3 



X: Exp't CENTER LINE VELOCITY, O: Exp't HALF JET WIDTH, ; 

COMPUT. CENTER LINE VELOCITY, COMPUT. HALF JET 

WIDTH 


Figure 6. A plane jet exhausting into a moving stream. 


11 


= 0.650 (X/d = 68) X = 0.381 (X/h = 40) 
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AND COMPUT. VELOCITY, TURBULENT KINETIC ENERGY, AND —u'v' RESPECTIVELY 




5. A Circular Jet Exhausting Into a Moving Stream 

Anomaly between plane jet and circular jet has been discussed more than enough 
number of times from the early days of k-e type turbulence model development. 
Recently, the anomaly has been somewhat resolved in Hanjelic and Launder [12] by 
introducing irrotational strain into the turbulence equations. Unfortunately, the 
improvement cannot be carried over to computation of multi- dimensional elliptic flows. 
The advantage of the present turbulence model also lies in the fact that the improve- 
ment can be carried over to elliptic flow computations, as can be found in Reference 7. 

A circular jet experiment due to Antonia and Bilger [23] is considered herein. 
The turbulent kinetic ener gy profile for initial condition data was prepared from the 

experimental data of /(u* )/Uq, where u’ is the fluctuating velocity in the flow direc- 
tion and Uq is the center-line excess velocity, by assuming that the turbulent kinetic 

energy would be approximately 1.2 times the flow direction normal Reynolds stress. 

The computational domain in the transverse direction extends from y = 0 m, the 
center-line, to y = 0.084 m, which corresponds to approximately nine times the half 
jet width. The flow direction domain starting from x = 0.2 m (x/d = 40, where d = 
0.00528 m is the diameter of the jet) to x = 1.4 m (x/d = 265) was discretized into 
1200 line-steps. An average of six iterations were required for each line-step to 
achieve the convergence criterion. 

Computational results showed that the flow field re-developed due to the 
inaccurate turbulent kinetic energy data; and that a similarity state which is slightly 
different from experimental data was achieved in the far down-stream region. Decay 
of the center-line velocity and the growth of the half jet width are shown in Figure 8, 
where it can be seen that the spreading rate at the far down- stream is close to that 
of the experimental data. 


cr* *5! 
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NOTATIONS ARE THE SAME AS IN FIGURE 2. 


Figure 8. A circular jet exhausting into a moving stream. 
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Due to the re-development of the flow field, the computational results are com- 
pared with experimental data in similarity coordinate (Fig. 9). It can be seen in 
Figure 9 that the velocity profile is flatter than that of the experimental data in the 
center region of the jet, and that the gradient of velocity profile is stiffer than 
experimental data in the middle region of the jet. The Reynolds stress profile 
obtained by using the Boussinesq eddy viscosity assumption that -u’v’ = v^9u/9y is 

higher than experimental data due to the stiffer velocity gradient in the same middle 
region. These discrepancies are due to the deficiency of the c^ function, and these 

degenerated computational results were not obtained in the finite difference computa- 
tion of the same flow case using a constant c^ [7]. 


6. A Wall Jet Issuing into a Moving Stream 

The wall jet is a boundary layer flow which finds a great number of applications 
in engineering processes. It is used in many different situations of film cooling pro- 
cesses and in preventing separation of boundary layer flows. The wall jet provides 
a serious test case for a numerical method as well as a turbulence model, since both 
of these must be able to predict the behavior of wall bounded boundary layer flows 
and jets separately. Finite difference computation of wall jet flows can be found in 
Ljuboja and Rodi [24] among many others. 

One of the most complete experimental data for the wall jet flows can be found 
in Irwin [25]. A configuration of wall jet is shown in Figure 10, where b is the jet 
slot height, Uj is the jet exit velocity, u^Cx) is the free stream velocity, Uq is the 

maximum excess velocity, u^(u = u„+u„) is the maximum velocity, y^ is a distance 

from the wall to a location where u = u^^^, and y-^j 2 3®^ width defined as 

the distance from the wall to the outer side location where u = UQ/2+Ug. Input data 

used in computation of the flow were obtained directly and/or by curve-fitting the 
experimental data [25]. These are: b = 0.00673 m, Uj = 60.64 m/sec; and the 

external free stream velocity and the length -scale of the flow field are given below: 

1 -fl 448 

Ue(x) = Uj [0.12052 (x/b) + 0.98425] (11) 

Yo^x) 

— = 0.04394 (x/b) + 0.3819 . (12) 


The computational domain extends from y = 0.002 m to y = 0.08 m in the trans- 
verse direction; and from x = 0.5532 m (x/b = 82.2) to x = 1.6892 m (x/b = 251) in 
the flow direction. The flow direction domain was discretized by 1135 line- steps, and 
an average of 14 iterations were required for each line- step to achieve the same 
convergence criterion given previously. It can be seen in Figures 11 and 12 that 
computational results compare favorably with experimental data in every detail. 
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Figure 9. A circular jet exhausting into a moving stream. 










CONCLUSIONS 


An algebraic stress turbulence model was presented together with finite element 
computation of a few turbulent flows. Experimentally deduced c^ curves reveal that 

the value of c^ needs to be close to that of the equilibrium condition (i.e., the pro- 
duction rate is the same as the dissipation rate) at the center region of jets where 
production rate is small but the turbulent kinetic energy is high. The center region 
of jets is apparently a high turbulence region even though the production rate is 
small. It is believed that further improvement of the c^ function to incorporate the 

experimentally observed behavior of c^ at the center region of jets will undoubtedly 

improve the turbulence model proposed. This fact can be further supported by the 
finite difference computation of the same turbulent boundary layer flows using the 
same turbulence model but with a constant c^ [7], where significantly improved com- 
putational results over the standard k-e turbulence model can be found not only for 
turbulent boundary layer flows but also for separated and swirling elliptic flows. In 
the present finite element computation, use of the c^ function was indispensable to 

achieve realizability (i.e., the turbulent kinetic energy and the dissipation rate cannot 
become negative). In fact, most of the k-e type turbulence models were able to run 
in the present boundary layer flow code [3] without any numerical difficulty if the 
c function was included into those turbulence models. In this sense, the present 

study clarified one of the fundamental differences between the finite element method 
and the finite difference method. 
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